/*

File name		: figure4.do

Description		: Grid search of the parameter γ
									
Output File		: figure4.png
				
*/

set more off, perm
capture restore

*** 001 Grid Search of parameter gamma ***
use "data/blood_donation/blood_donation.dta", clear

	// Merge rainfall data
	merge m:1 entnahmedatum using "data/rainfall/rainfallZH.dta"
	*, keep(match)

	// Rename kalendertag daily_rainfall
	rename daily_rainfall_ZH daily_rainfall
	gen any_rain = 1 if daily_rainfall > 0 & daily_rainfall != .
	replace any_rain = 0 if daily_rainfall == 0
	sort spnr entnahmedatum
	gen any_rain_plus1 = any_rain[_n+1] if spnr == spnr[_n+1]

	sort spnr period
	sum don if period == 1
	sum don if period == 2 & don[_n-1] == 1 & spnr == spnr[_n-1]
	sum don if period == 2 & don[_n-1] == 0 & spnr == spnr[_n-1]
	sum don if period == 3 & don[_n-1] == 1 & spnr == spnr[_n-1]
	sum don if period == 3 & don[_n-1] == 0 & spnr == spnr[_n-1]
	sum don if period == 4 & don[_n-1] == 1 & spnr == spnr[_n-1]
	sum don if period == 4 & don[_n-1] == 0 & spnr == spnr[_n-1]

	/* Gen new variable that is 1 = donated in period 1, 2 = did not donate in period 1, 
	 3 = donated in period2, 4 = did not donate in period 2 */
	sort spnr period
	gen treatment_bydonation1 = 1 if period == 1 & don == 1
	replace treatment_bydonation1 = 0 if period == 1 & don == 0
	replace treatment_bydonation1 = treatment_bydonation1[_n-1] if spnr == spnr[_n-1]
	gen treatment_bydonation2 = 1 if period == 2 & don == 1
	replace treatment_bydonation2 = 0 if period == 2 & don == 0
	replace treatment_bydonation2 = treatment_bydonation2[_n-1] ///
	if spnr == spnr[_n-1] & treatment_bydonation2[_n-1]!= .
	replace treatment_bydonation2 = treatment_bydonation2[_n+1] ///
	if spnr == spnr[_n+1] & treatment_bydonation2[_n+1]!= .

// Define grid search program
global within_type .
capture program drop within2
program within2
	
	local idfe_yes $idfe_yes
	
	capture drop loglik1 loglik2
	capture drop loglikk1 loglikk2
	capture drop reached_hat reached_minus1_hat reached_minus2_hat reached_minus3_hat
	
	matrix ll = (.,.)
	matrix lll = (.,.)

	foreach l of varlist reached reached_minus1 reached_minus2 reached_minus3 {
		
		if `idfe_yes' == 0 {
				qui ivreghdfe `l' phonecall phonecall_minus1 phonecall_minus2 ///
				phonecall_minus3 male age bltype_dummy2 bltype_dummy3, ///
				absorb(c_num)
			}
			else{
				qui ivreghdfe `l' phonecall phonecall_minus1 phonecall_minus2 ///
				phonecall_minus3, ///
				absorb(c_num spnr)
			}
		
		predict double `l'_hat, xb
		
	}


	forvalues i = 0.1(0.01)0.9 {
		capture drop reached_vector
		gen reached_vector = reached_hat + `i' * reached_minus1_hat ///
		+ `i'^2 * reached_minus2_hat + `i'^3 * reached_minus3_hat
		
		if `idfe_yes' == 0 {
				qui ivreghdfe don reached_vector ///
				male age bltype_dummy2 bltype_dummy3, ///
				absorb(c_num)
			}
			else{
				qui ivreghdfe don reached_vector, ///
				absorb(c_num spnr)
			}
		
			matrix ll[1,1] = `i'
			matrix ll[1,2] = `e(rss)'
		
		if `i' == 0.1 {
			matrix loglik = ll
		}
		else{
			matrix loglik = (loglik \ ll) 
		}
	}
	
	*max
	svmat loglik
	sort loglik2 loglik1
	sum loglik1 if _n == 1 & loglik1[1]>0.1 & loglik1[1]<0.9
	
	local max = `r(mean)'
	local lb = `max'-0.01
	local ub = `max'+0.01
	
	forvalues j = `lb'(0.001)`ub' {
		
		capture drop reached_vector
		gen reached_vector = reached_hat + `j' * reached_minus1_hat + ///
		`j'^2 * reached_minus2_hat + `j'^3 * reached_minus3_hat
		
		if `idfe_yes' == 0 {
				qui ivreghdfe don reached_vector ///
				male age bltype_dummy2 bltype_dummy3, ///
				absorb(c_num)
			}
			else{
				qui ivreghdfe don reached_vector, ///
				absorb(c_num spnr)
			}
			
			matrix lll[1,1] = `j'
			matrix lll[1,2] = `e(rss)'
		
		if `j' == `lb'{
			matrix loglikk = lll
		}
		else{
			matrix loglikk = (loglikk \ lll) 
		}

	}
	
	*max
	svmat loglikk
	sort loglikk2 loglikk1
	matrix max = .
	sum loglikk1 if _n == 1 & loglikk1[1] > 0.1 & loglikk1[1] < 0.9
	local rounded = round(`r(mean)', .001)
	capture matrix max = `rounded'

end

// Run within2 first without individual fixed effects
global idfe_yes 0
	within2
	matrix maxi = max
	gen loglik2_fig = loglik2
	gen loglik1_fig = loglik1

// Then run within2 first with individual fixed effects
global idfe_yes 1
	within2

	graph twoway (scatter loglik2_fig loglik1_fig, color(navy) yaxis(2) ///
	ylabel(758.5 (0.5)760.5, axis(2)) yscale(axis(2) range(758.5 760.5) ///
	titlegap(3)) text(427.77 0.1 "Blood drive FE" "& individual FE", ///
	place(n) color(navy)) text(427.43 0.1 "Blood drive FE", place(n) color(navy)) ///
	text(427.77 0.484 "{&gamma} = 0.484"" (0.131)", place(n) color(maroon)) ///
	text(427.89 0.484 "{&beta}{subscript:1} = 0.175"" (0.030)", place(n) color(black))) ///
	(scatter loglik2 loglik1, yaxis(1) color(navy) ///
	graphregion(color(white)) xscale(titlegap(3)) ylabel(426.7 (0.3) 427.9, axis(1)) ///
	yscale(axis(1) range(426.7 427.9) titlegap(3)) bgcolor(white) ///
	ytitle("Residual sum of squares", axis(1)) ytitle("Residual sum of squares", axis(2)) ///
	xtitle("Habit formation parameter") legend(off) ///
	text(427.77 0.692 "{&gamma} = 0.692"" (0.189)", place(n) color(maroon)) ///
	text(427.89 0.692 "{&beta}{subscript:1} = 0.246"" (0.045)", place(n) color(black))) ///
	(scatteri 426.7 0.692 427.75 0.692, color(maroon) yaxis(1) recast(line) lpattern(dash)) ///
	(scatteri 758.5 0.484 760.25 0.484, color(maroon) yaxis(2) recast(line) lpattern(dash))
	
	graph export output/figures/figure4.png, replace width(2048)	
